#include"mrrr.h"
/* Subroutine */int pdlasq4(int *i0, int *n0, double *z__, int *pp, int *n0in,
		double *dmin__, double *dmin1, double *dmin2, double *dn, double *dn1,
		double *dn2, double *tau, int *ttype, double *g) {
	/* System generated locals */
	int i__1;
	double d__1, d__2;


	/* Local variables */
	double s, a2, b1, b2;
	int i4, nn, np;
	double gam, gap1, gap2;

	/*  -- LAPACK routine (version 3.2)                                    -- */

	/*  -- Contributed by Osni Marques of the Lawrence Berkeley National   -- */
	/*  -- Laboratory and Beresford Parlett of the Univ. of California at  -- */
	/*  -- Berkeley                                                        -- */
	/*  -- November 2008                                                   -- */

	/*  -- LAPACK is a software package provided by Univ. of Tennessee,    -- */
	/*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- */

	/*     .. Scalar Arguments .. */
	/*     .. */
	/*     .. Array Arguments .. */
	/*     .. */

	/*  Purpose */
	/*  ======= */

	/*  DLASQ4 computes an approximation TAU to the smallest eigenvalue */
	/*  using values of d from the previous transform. */

	/*  I0    (input) int */
	/*        First index. */

	/*  N0    (input) int */
	/*        Last index. */

	/*  Z     (input) DOUBLE PRECISION array, dimension ( 4*N ) */
	/*        Z holds the qd array. */

	/*  PP    (input) int */
	/*        PP=0 for ping, PP=1 for pong. */

	/*  NOIN  (input) int */
	/*        The value of N0 at start of EIGTEST. */

	/*  DMIN  (input) DOUBLE PRECISION */
	/*        Minimum value of d. */

	/*  DMIN1 (input) DOUBLE PRECISION */
	/*        Minimum value of d, excluding D( N0 ). */

	/*  DMIN2 (input) DOUBLE PRECISION */
	/*        Minimum value of d, excluding D( N0 ) and D( N0-1 ). */

	/*  DN    (input) DOUBLE PRECISION */
	/*        d(N) */

	/*  DN1   (input) DOUBLE PRECISION */
	/*        d(N-1) */

	/*  DN2   (input) DOUBLE PRECISION */
	/*        d(N-2) */

	/*  TAU   (output) DOUBLE PRECISION */
	/*        This is the shift. */

	/*  TTYPE (output) int */
	/*        Shift type. */

	/*  G     (input/output) REAL */
	/*        G is passed as an argument in order to save its value between */
	/*        calls to DLASQ4. */

	/*  Further Details */
	/*  =============== */
	/*  CNST1 = 9/16 */

	/*  ===================================================================== */

	/*     .. Parameters .. */
	/*     .. */
	/*     .. Local Scalars .. */
	/*     .. */
	/*     .. Intrinsic Functions .. */
	/*     .. */
	/*     .. Executable Statements .. */

	/*     A negative DMIN forces the shift to take that absolute value */
	/*     TTYPE records the type of shift. */

	/* Parameter adjustments */
	--z__;

	/* Function Body */
	if (*dmin__ <= 0.) {
		*tau = -(*dmin__);
		*ttype = -1;
		return 0;
	}

	nn = (*n0 << 2) + *pp;
	if (*n0in == *n0) {

		/*        No eigenvalues deflated. */

		if (*dmin__ == *dn || *dmin__ == *dn1) {

			b1 = sqrt(z__[nn - 3]) * sqrt(z__[nn - 5]);
			b2 = sqrt(z__[nn - 7]) * sqrt(z__[nn - 9]);
			a2 = z__[nn - 7] + z__[nn - 5];

			/*           Cases 2 and 3. */

			if (*dmin__ == *dn && *dmin1 == *dn1) {
				gap2 = *dmin2 - a2 - *dmin2 * .25;
				if (gap2 > 0. && gap2 > b2) {
					gap1 = a2 - *dn - b2 / gap2 * b2;
				} else {
					gap1 = a2 - *dn - (b1 + b2);
				}
				if (gap1 > 0. && gap1 > b1) {
					/* Computing MAX */
					d__1 = *dn - b1 / gap1 * b1, d__2 = *dmin__ * .5;
					s = max(d__1, d__2);
					*ttype = -2;
				} else {
					s = 0.;
					if (*dn > b1) {
						s = *dn - b1;
					}
					if (a2 > b1 + b2) {
						/* Computing MIN */
						d__1 = s, d__2 = a2 - (b1 + b2);
						s = min(d__1, d__2);
					}
					/* Computing MAX */
					d__1 = s, d__2 = *dmin__ * .333;
					s = max(d__1, d__2);
					*ttype = -3;
				}
			} else {

				/*              Case 4. */

				*ttype = -4;
				s = *dmin__ * .25;
				if (*dmin__ == *dn) {
					gam = *dn;
					a2 = 0.;
					if (z__[nn - 5] > z__[nn - 7]) {
						return 0;
					}
					b2 = z__[nn - 5] / z__[nn - 7];
					np = nn - 9;
				} else {
					np = nn - (*pp << 1);
					b2 = z__[np - 2];
					gam = *dn1;
					if (z__[np - 4] > z__[np - 2]) {
						return 0;
					}
					a2 = z__[np - 4] / z__[np - 2];
					if (z__[nn - 9] > z__[nn - 11]) {
						return 0;
					}
					b2 = z__[nn - 9] / z__[nn - 11];
					np = nn - 13;
				}

				/*              Approximate contribution to norm squared from I < NN-1. */

				a2 += b2;
				i__1 = (*i0 << 2) - 1 + *pp;
				for (i4 = np; i4 >= i__1; i4 += -4) {
					if (b2 == 0.) {
						goto L20;
					}
					b1 = b2;
					if (z__[i4] > z__[i4 - 2]) {
						return 0;
					}
					b2 *= z__[i4] / z__[i4 - 2];
					a2 += b2;
					if (max(b2, b1) * 100. < a2 || .563 < a2) {
						goto L20;
					}
					/* L10: */
				}
				L20: a2 *= 1.05;

				/*              Rayleigh quotient residual bound. */

				if (a2 < .563) {
					s = gam * (1. - sqrt(a2)) / (a2 + 1.);
				}
			}
		} else if (*dmin__ == *dn2) {

			/*           Case 5. */

			*ttype = -5;
			s = *dmin__ * .25;

			/*           Compute contribution to norm squared from I > NN-2. */

			np = nn - (*pp << 1);
			b1 = z__[np - 2];
			b2 = z__[np - 6];
			gam = *dn2;
			if (z__[np - 8] > b2 || z__[np - 4] > b1) {
				return 0;
			}
			a2 = z__[np - 8] / b2 * (z__[np - 4] / b1 + 1.);

			/*           Approximate contribution to norm squared from I < NN-2. */

			if (*n0 - *i0 > 2) {
				b2 = z__[nn - 13] / z__[nn - 15];
				a2 += b2;
				i__1 = (*i0 << 2) - 1 + *pp;
				for (i4 = nn - 17; i4 >= i__1; i4 += -4) {
					if (b2 == 0.) {
						goto L40;
					}
					b1 = b2;
					if (z__[i4] > z__[i4 - 2]) {
						return 0;
					}
					b2 *= z__[i4] / z__[i4 - 2];
					a2 += b2;
					if (max(b2, b1) * 100. < a2 || .563 < a2) {
						goto L40;
					}
					/* L30: */
				}
				L40: a2 *= 1.05;
			}

			if (a2 < .563) {
				s = gam * (1. - sqrt(a2)) / (a2 + 1.);
			}
		} else {

			/*           Case 6, no information to guide us. */

			if (*ttype == -6) {
				*g += (1. - *g) * .333;
			} else if (*ttype == -18) {
				*g = .083250000000000005;
			} else {
				*g = .25;
			}
			s = *g * *dmin__;
			*ttype = -6;
		}

	} else if (*n0in == *n0 + 1) {

		/*        One eigenvalue just deflated. Use DMIN1, DN1 for DMIN and DN. */

		if (*dmin1 == *dn1 && *dmin2 == *dn2) {

			/*           Cases 7 and 8. */

			*ttype = -7;
			s = *dmin1 * .333;
			if (z__[nn - 5] > z__[nn - 7]) {
				return 0;
			}
			b1 = z__[nn - 5] / z__[nn - 7];
			b2 = b1;
			if (b2 == 0.) {
				goto L60;
			}
			i__1 = (*i0 << 2) - 1 + *pp;
			for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
				a2 = b1;
				if (z__[i4] > z__[i4 - 2]) {
					return 0;
				}
				b1 *= z__[i4] / z__[i4 - 2];
				b2 += b1;
				if (max(b1, a2) * 100. < b2) {
					goto L60;
				}
				/* L50: */
			}
			L60: b2 = sqrt(b2 * 1.05);
			/* Computing 2nd power */
			d__1 = b2;
			a2 = *dmin1 / (d__1 * d__1 + 1.);
			gap2 = *dmin2 * .5 - a2;
			if (gap2 > 0. && gap2 > b2 * a2) {
				/* Computing MAX */
				d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
				s = max(d__1, d__2);
			} else {
				/* Computing MAX */
				d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
				s = max(d__1, d__2);
				*ttype = -8;
			}
		} else {

			/*           Case 9. */

			s = *dmin1 * .25;
			if (*dmin1 == *dn1) {
				s = *dmin1 * .5;
			}
			*ttype = -9;
		}

	} else if (*n0in == *n0 + 2) {

		/*        Two eigenvalues deflated. Use DMIN2, DN2 for DMIN and DN. */

		/*        Cases 10 and 11. */

		if (*dmin2 == *dn2 && z__[nn - 5] * 2. < z__[nn - 7]) {
			*ttype = -10;
			s = *dmin2 * .333;
			if (z__[nn - 5] > z__[nn - 7]) {
				return 0;
			}
			b1 = z__[nn - 5] / z__[nn - 7];
			b2 = b1;
			if (b2 == 0.) {
				goto L80;
			}
			i__1 = (*i0 << 2) - 1 + *pp;
			for (i4 = (*n0 << 2) - 9 + *pp; i4 >= i__1; i4 += -4) {
				if (z__[i4] > z__[i4 - 2]) {
					return 0;
				}
				b1 *= z__[i4] / z__[i4 - 2];
				b2 += b1;
				if (b1 * 100. < b2) {
					goto L80;
				}
				/* L70: */
			}
			L80: b2 = sqrt(b2 * 1.05);
			/* Computing 2nd power */
			d__1 = b2;
			a2 = *dmin2 / (d__1 * d__1 + 1.);
			gap2 = z__[nn - 7] + z__[nn - 9]
					- sqrt(z__[nn - 11]) * sqrt(z__[nn - 9]) - a2;
			if (gap2 > 0. && gap2 > b2 * a2) {
				/* Computing MAX */
				d__1 = s, d__2 = a2 * (1. - a2 * 1.01 * (b2 / gap2) * b2);
				s = max(d__1, d__2);
			} else {
				/* Computing MAX */
				d__1 = s, d__2 = a2 * (1. - b2 * 1.01);
				s = max(d__1, d__2);
			}
		} else {
			s = *dmin2 * .25;
			*ttype = -11;
		}
	} else if (*n0in > *n0 + 2) {

		/*        Case 12, more than two eigenvalues deflated. No information. */

		s = 0.;
		*ttype = -12;
	}

	*tau = s;
	return 0;

	/*     End of DLASQ4 */

} /* dlasq4_ */
